#!/usr/local/bin/python

##########################################################################
#
#  pom2k.py -- Python port of the Princeton Ocean Model.
#  Copyright (C) 2013 Konstantin Grigoryev
#  
#  This program is free software; you can redistribute it and/or modify
#  it under the terms of the GNU General Public License as published by
#  the Free Software Foundation; either version 2 of the License, or
#  (at your option) any later version.
#  
#  This program is distributed in the hope that it will be useful,
#  but WITHOUT ANY WARRANTY; without even the implied warranty of
#  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
#  GNU General Public License for more details.
#  
#  You should have received a copy of the GNU General Public License along
#  with this program; if not, write to the Free Software Foundation, Inc.,
#  51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
#
##########################################################################


# This code is based on Fortran POM sources, last changed on 2006-05-03,
# which is yet available at http://www.ccpo.odu.edu/POMWEB/POM2K/
# Authors of the original code: A. Blumberg and G. Mellor.
# Author of this Python port: K. Grigoryev.
#
# The following note was copied from the original code, excluding some
# Fortran-specific sections.
#
#
# FUNCTION
#
# This is a version of the three dimensional, time dependent, primitive
# equation, ocean model developed by Alan Blumberg and George Mellor with
# subsequent contributions by Leo Oey, Steve Brenner and others. It is now
# called the Princeton Ocean Model. Two references are:
# 
#   Blumberg, A.F. and G.L. Mellor; Diagnostic and prognostic numerical
#    circulation studies of the South Atlantic Bight, J. Geophys. Res. 88,
#    4579-4592, 1983.
#
#   Blumberg, A.F. and G.L. Mellor; A description of a three-dimensional
#    coastal ocean circulation model, Three-Dimensional Coastal Ocean
#    Models, Coastal and Estuarine Sciences, 4, N.S. Heaps, ed., American
#    Geophysical Union, 1-16, 1987.
#
# In subroutine profq the model makes use of the turbulence closure
# sub-model described in:
#
#   Mellor, G.L. and T. Yamada; Development of a turbulence closure model
#    for geophysical fluid problems, Rev. Geophys. Space Phys., 20, No. 4,
#    851-875, 1982.
# (note recent profq that includes breaking waves)
#
# A user's guide is available:
#
#   Mellor, G.L.; User's guide for a three-dimensional, primitive
#    equation, numerical ocean model. Princeton University Report, 1998.
#
# In October 2001, the source code underwent revision by John Hunter of
# the University of Tasmania. Major aspects of the revision were:
#
#  (1) The revision was based on pom98 updated to 12/9/2001.
#  (2) Declaration of all variables.
#  (3) Rationalisation of the input of all constants.
#  (4) Modifications to the "printer" output.
#  (5) Output to a netCDF file.
#  (6) Inclusion of surface freshwater flux.
#  (7) Inclusion of atmospheric pressure.
#  (8) Inclusion of an additional problem to check (6) and (7), above.
#  (9) Inclusion of option for Smolarkiewicz advection scheme.
#
# This revised version is functionally almost equivalent to pom98. The
# output to device 6 from the "seamount" problem should be almost the
# same, any differences being due to minor format changes and improvements
# in rounding.
#
# This revision was helped by the following people: Tal Ezer, Peter
# Holloway, George Mellor, Rich Signell, Ian Webster, Brian Williams
# and Emma Young.
#
#
# GENERAL NOTES
#
# 1. All units are S.I. (M.K.S.) unless otherwise stated. NOTE that time
#    is in days from the start of the run.
#
# 2. "b", <nothing> and "f" refers to backward, central and forward time
#    levels.
#
# 3. NetCDF output may be used. In order to omit/use netCDF, comment/
#    /uncomment all statements carrying the comment "*netCDF*" at the end
#    of the line (or set netcdf_file='nonetcdf')
#
# 4. NetCDF is version 3. An attempt has been made to conform to the
#    NetCDF Climate and Forecast (CF) Metadata Conventions, but this may
#    not yet be complete.
#    (see: http://www.cgd.ucar.edu/cms/eaton/cf-metadata/index.html)
#
# 5. In order to use netCDF, the [original -- K.G.] program should be
#    compiled with the appropriate library. For example, if using g77, you
#    may need to type:
#
#      g77 -o pom2k pom2k.f /usr/lib/libnetcdf.a
# 
#    You should also have the "include" file of netCDF subroutines
#    (pom2k.n).
#
# 6. In order to use netCDF, you may need to change the name of the
#    "include" file in the statement:
#
#      include '/usr/include/netcdf.inc'
#
#    in subroutine write_netcdf.


from numpy import *

# In the original program, 'source' should be agreed with 'source_c'
# in pom2k.c and 'source_n' in pom2k.n. Here it's just a version mark.
source = 'pom2k  2006-05-03'
 
# Set size of problem here.
(im, jm, kb) = (41, 61, 16)     ##TODO Add loading from the 'grid' file

small = 1.e-9           # Small value, used to avoid division by zero
spdy  = 86400.0         # Day duration, in seconds ("spdy" means "secs-per-day")
#pi    = 4 * atan(1)    # Pi. Commented out, because it's imported from 'numpy'

title = 'Run 1'         # Run's title, for pretty output only

# Calculation mode:
#   1. Not allowed;
#   2. 2-D calculation (bottom stress calculated in advave);
#   3. 3-D calculation (bottom stress calculated in profu,v);
#   4. 3-D calculation with t and s held fixed.
mode = 3                ##TODO Magic number!

# Advection scheme:
#   1. Centred scheme, as originally provide in POM;
#   2. Smolarkiewicz iterative upstream scheme, based on subroutines
#      provided by Gianmaria Sannino and Vincenzo Artale.
nadv = 1                ##TODO Magic number!

# Constants for Smolarkiewicz iterative upstream scheme. Number of
# iterations. This should be in the range 1 - 4. 1 is standard upstream
# differencing; 3 adds 50% CPU time to POM.
nitera = 2

# Smoothing parameter. This should preferably be 1, but 0 < sw < 1 gives
# smoother solutions with less overshoot when nitera > 1.
sw = 0.5

# Index to indicate whether run to start from restart file.
#   0. No restart input file;
#   1. Restart input file.
nread = 0               ##TODO Magic number!

# External (2-D) time step (secs.) according to CFL.
dte = 6.0

# Internal (3-D) time step>/<External (2-D) time step>
# (dti/dte_) -- dimensionless, integer.
isplit = 30

# Date and time of start of initial run of model in format (i.e.
# UDUNITS convention)
#
#   YYYY-MM-DD HH:MM:SS <+/->HH:MM
#
# where "<+/->HH:MM" is the time zone (positive eastwards from
# Coordinated Universal Time). Note that the climatological time axis
# (i.e. beginning of year zero, which does not exist in the real-world
# calendar) has been used here. Insert your own date and time as required.
time_start = '2000-01-01 00:00:00 +00:00'

days  = 0.25            # Run duration in days
prtd1 = 0.125           # Initial print interval in days
prtd2 = 1.0             # Final print interval in days
swtch = 1000.0          # Time to switch from prtd1 to prtd2
iskp  = 4               # Printout skip interval in i
jskp  = 3               # Printout skip interval in j

# Logical for inertial ramp (True if inertial ramp to be applied to wind
# stress and baroclinic forcing, otherwise False)
lramp = False           ##TODO Magic number!

# Reference density. Recommended values: 1025 for seawater, 1000 for
# freswater -- S.I. units.
rhoref = 1025.0

tbias  = 0.0            # Temperature bias -- deg.C
sbias  = 0.0            # Salinity bias
grav   = 9.80665        # Gravity constant -- S.I. units
kappa  = 0.4            # Von Karman's constant
z0b    = 0.01           # Bottom roughness -- metres
cbcmin = 0.0025         # Minimum bottom friction coef.
cbcmax = 1.0            # Maximum bottom friction coef.
horcon = 0.2            # Smagorinsky diffusivity coef.

# Inverse horizontal turbulent Prandtl number (ah/am) -- dimensionless.
# Nhat tprni = 0 yields zero horizontal diffusivity!
tprni = 0.2

# Background viscosity used in subroutines profq, proft, profu and
# profv -- S.I. units.
umol = 2.e-5

# Maximum depth used in radiation boundary condition in subroutine
# bcond -- metres.
hmax = 4500.0

# Maximum magnitude of vaf (used in check that essentially tests for CFL
# violation).
vmaxl = 100.0

# Maximum allowable value of <difference of depths>/<sum of depths> for
# two adjacent cells (dimensionless). This is used in subroutine slpmax.
# If >= 1, then slpmax is not applied.
slmax = 2.0

# Jerlov water type:
#   ntp :   1    2    3    4    5
#   type:   I  I-a  I-b   II  III
ntp = 2                 ##TODO Magic number!

# Surface temperature boundary condition, used in subroutine proft:
#   1. Prescribed flux, short wave penetration: no;
#   2. Prescribed flux, short wave penetration: yes;
#   3. Prescribed temperature, short wave penetration: no;
#   4. Prescribed temperature, short wave penetration: yes.
nbct = 1                ##TODO Magic number!

# Surface salinity boundary condition, used in subroutine proft:
#   1. Prescribed flux;
#   2. Not allowed;
#   3. Prescribed salinity;
#   4. Not allowed.
nbcs = 1                ##TODO Magic number!

# Step interval during which external (2-D) mode advective terms are
# not updated -- dimensionless.
ispadv = 5

# Constant in temporal filter used to prevent solution splitting --
# dimensionless.
smoth = 0.1

# Weight used for surface slope term in external (2-D) dynamic
# equation (a value of alpha = 0 is perfectly acceptable, but the
# value, alpha = 0.225 permits a longer time step).
alpha = 0.225

# Initial value of aam.
aam_init = 500.0


# Sigma coordinate which spans the domain, z = 0 (surface) to
# z = -1 (bottom).
z = empty(kb)

# Sigma coordinate, intermediate between z.
zz = empty(kb)

# dz(k) = z(k) - z(k+1),  dzz(k) = zz(k) - zz(k+1)
dz  = empty(kb)
dzz = empty(kb)

# Vertical average of aam -- m2 / s.
aam2d = empty((im, jm))

# Sum of the second, third and fourth terms in equations (18,19).
advua = empty((im, jm))
advva = empty((im, jm))

# Vertical integrals of advx, advy; also the sum of the fourth, fifth and
# sixth terms in equations (22a,b).
adx2d = empty((im, jm))  
ady2d = empty((im, jm))

# Cell areas centered on the variables, T, U and V respectively -- m2.
art = empty((im, jm))
aru = empty((im, jm))
arv = empty((im, jm))

# Bottom friction coefficient.
cbc = empty((im, jm))

# The Coriolis parameter -- 1/s.
cor = empty((im, jm))

# d = h + el,  dt = h + et.
d  = empty((im, jm))
dt = empty((im, jm))

# Vertical integrals of drhox and drhoy.
drx2d = empty((im, jm))  
dry2d = empty((im, jm))

# Mask for scalar variables; 0 over land, 1 over water.
fsm = empty((im, jm), dtype=byte)

# Masks for the u and v component of velocity; 0 over land, 1 over water.
dum = empty((im, jm), dtype=byte)
dvm = empty((im, jm), dtype=byte)

# Grid spacing -- metres.
dx = empty((im, jm))
dy = empty((im, jm))

east_c  = empty((im, jm))
east_e  = empty((im, jm))
east_u  = empty((im, jm))
east_v  = empty((im, jm))

north_c = empty((im, jm))
north_e = empty((im, jm))
north_u = empty((im, jm))
north_v = empty((im, jm))

# Atmospheric pressure.
e_atmos = empty((im, jm))

# The surface elevation asused in the internal mode for the pressure
# gradient and derived from el -- metres.
egb = empty((im, jm))    
egf = empty((im, jm))

# The surface elevation as used in the external mode -- metres.
el  = empty((im, jm))
elb = empty((im, jm))
elf = empty((im, jm))

# The surface elevation as used in the internal mode and derived
# from el -- metres.
et  = empty((im, jm))
etb = empty((im, jm))
etf = empty((im, jm))

fluxua = empty((im, jm)) 
fluxva = empty((im, jm))

# The bottom depth -- metres.
h = empty((im, jm))

psi = empty((im, jm))
rot = empty((im, jm))

# Temperature and salinity at the surface (?) -- deg.C, psu.
tsurf = empty((im, jm))
ssurf = empty((im, jm))

# Short wave radiation incident on the ocean surface -- deg.C m / s.
swrad = empty((im, jm))

vfluxb = empty((im, jm)) 
vfluxf = empty((im, jm))

# External (2-D) CFL time step.
tps = empty((im, jm))

# Vertical mean of U, V -- metres per second.
ua  = empty((im, jm))
uab = empty((im, jm))
uaf = empty((im, jm))
va  = empty((im, jm))
vab = empty((im, jm))
vaf = empty((im, jm))

# ua, va time averaged over the interval -- metres per second.
utb = empty((im, jm))
utf = empty((im, jm))
vtb = empty((im, jm))
vtf = empty((im, jm))

# Temperature and salinity fluxes at the surface -- ged.C m / s, psu m / s.
wtsurf = empty((im, jm))
wssurf = empty((im, jm))

# Momentum fluxes at the surface and at the bottom -- m2 / s2.
wusurf = empty((im, jm))
wvsurf = empty((im, jm))
wubot  = empty((im, jm))
wvbot  = empty((im, jm))

# Horizontal kinematic viscosity -- m2 / s.
aam = empty((im,jm,kb))

# Horizontal advection and diffusion terms in equations (3) and (4).
advx = empty((im,jm,kb))
advy = empty((im,jm,kb))

a = empty((im,jm,kb))
c = empty((im,jm,kb))

# Internal baroclinic pressure gradient, x and y components.
drhox = empty((im,jm,kb))
drhoy = empty((im,jm,kb))

dtef = empty((im,jm,kb)) 

ee = empty((im,jm,kb))
gg = empty((im,jm,kb))

# Vertical diffusivity -- m2 / s.
kh = empty((im,jm,kb))

# Vertical kinematic viscosity -- m2 / s.
km = empty((im,jm,kb))

kq = empty((im,jm,kb))

# Turbulence length scale -- metres;
l = empty((im,jm,kb))

# Twice the turbulence kinetic energy -- m2 / s2.
q2  = empty((im,jm,kb))
q2b = empty((im,jm,kb))

# q2l = q2 * l.
q2l  = empty((im,jm,kb))
q2lb = empty((im,jm,kb))

# Relaive density -- dimensionless.
rho = empty((im,jm,kb))

# Density field which is horizontally averaged before transfer to sigma
# coordinates.
rmean = empty((im,jm,kb))

# Salinity -- psu.
s  = empty((im,jm,kb))
sb = empty((im,jm,kb))

# Temperature -- deg.C.
t  = empty((im,jm,kb))
tb = empty((im,jm,kb))

# Stationary silinity and temperature fields which approximately have
# the same vertical structure as s and t.
sclim = empty((im,jm,kb))
tclim = empty((im,jm,kb))

# Horizontal velocities -- metres per second.
u  = empty((im,jm,kb))
ub = empty((im,jm,kb))
uf = empty((im,jm,kb))
v  = empty((im,jm,kb))
vb = empty((im,jm,kb))
vf = empty((im,jm,kb))

# Sigma coordinate vertical velocity -- metres per second.
w = empty((im,jm,kb))

# Total vertical flux --
zflux = empty((im,jm,kb))



def daysToInternalSteps(days):
    return int(days * spdy / dti + 0.5)

dti  = dte * isplit
dte2 = 2 * dte          ##TODO Scope pollution            
dti2 = 2 * dti          ##TODO Scope pollution

iend   = max(2, daysToInternalSteps(days))
iprint = daysToInternalSteps(prtd1)
iswtch = daysToInternalSteps(swtch)

ispi  = 1.0 / isplit    ##TODO Scope pollution
isp2i = 0.5 / isplit    ##TODO Scope pollution

# Print initial summary
print
print ' source   = {0:<40}'.format(source)
print ' title      = {0:<40}'.format(title)
print ' iproblem   = {0:10d}'.format(3)     ##TODO Is this necessary?
print ' mode       = {0:10d}'.format(mode)
print ' nadv       = {0:10d}'.format(nadv)
print ' nitera     = {0:10d}'.format(nitera)
print ' sw         = {0:10.4f}'.format(sw)
print ' nread      = {0:10d}'.format(nread)   
print ' dte        = {0:10.2f}'.format(dte)
print ' dti        = {0:10.1f}'.format(dti)
print ' isplit     = {0:10d}'.format(isplit) 
print ' time_start = {0:<25}'.format(time_start)
print ' days       = {0:10.4f}'.format(days)
print ' iend       = {0:10d}'.format(iend)
print ' prtd1      = {0:10.4f}'.format(prtd1)
print ' iprint     = {0:10d}'.format(iprint)
print ' prtd2      = {0:10.4f}'.format(prtd2)
print ' swtch      = {0:10.2f}'.format(swtch)
print ' iswtch     = {0:10d}'.format(iswtch)
print ' iskp, jskp = {0:5d},{1:5d}'.format(iskp, jskp)
print ' lramp      = {0:>10s}'.format('T' if lramp else 'F')
print ' rhoref     = {0:10.3f}'.format(rhoref)
print ' tbias      = {0:10.3f}'.format(tbias)
print ' sbias      = {0:10.3f}'.format(sbias)
print ' grav       = {0:10.4f}'.format(grav)
print ' kappa      = {0:10.4f}'.format(kappa)
print ' z0b        = {0:10.6f}'.format(z0b)
print ' cbcmin     = {0:10.6f}'.format(cbcmin)
print ' cbcmax     = {0:10.6f}'.format(cbcmax)
print ' horcon     = {0:10.3f}'.format(horcon)
print ' tprni      = {0:10.4f}'.format(tprni)
print ' umol       = {0:10.4f}'.format(umol)
print ' hmax       = {0:10.2f}'.format(hmax)
print ' vmaxl      = {0:10.4f}'.format(vmaxl)
print ' slmax      = {0:10.4f}'.format(slmax)
print ' kl1, kl2   = {0:5d},{1:5d}'.format(kl1, kl2)
print ' ntp        = {0:10d}'.format(ntp)
print ' nbct       = {0:10d}'.format(nbct)
print ' nbcs       = {0:10d}'.format(nbcs)
print ' ispadv     = {0:10d}'.format(ispadv)
print ' smoth      = {0:10.4f}'.format(smoth)
print ' alpha      = {0:10.4f}'.format(alpha)

# Initialise boundary arrays
vabn[:] = 0.0
vabs[:] = 0.0
eln [:] = 0.0
els [:] = 0.0
vbn [:] = 0.0
vbs [:] = 0.0
tbn [:] = 0.0
tbs [:] = 0.0
sbn [:] = 0.0
sbs [:] = 0.0
uabe[:] = 0.0
uabw[:] = 0.0
ele [:] = 0.0
elw [:] = 0.0
ube [:] = 0.0
ubw [:] = 0.0
tbe [:] = 0.0
tbw [:] = 0.0
sbe [:] = 0.0
sbw [:] = 0.0

# Initialise 2-D and 3-D arrays for safety (this may be overwritten later)
uab    [:] = 0.0
vab    [:] = 0.0
elb    [:] = 0.0
etb    [:] = 0.0
e_atmos[:] = 0.0
vfluxb [:] = 0.0
vfluxf [:] = 0.0
wusurf [:] = 0.0
wvsurf [:] = 0.0
wtsurf [:] = 0.0
wssurf [:] = 0.0
swrad  [:] = 0.0
drx2d  [:] = 0.0
dry2d  [:] = 0.0
ub     [:] = 0.0
vb     [:] = 0.0

##
##TODO Here is a perfect place to define subroutines
##
    
# Read in grid data, and initial and lateral boundary conditions
file2ic()
 
# Inertial period for temporal filter
##FIXME (im/2 - 1, jm/2 - 1) is not a central point
period = (2 * pi) / abs(cor[im / 2 - 1, jm / 2 - 1]) / spdy

# Initialise time
time = time0 = 0.0

# Initial conditions. Note that lateral thermodynamic boundary conditions
# are often set equal to the initial conditions and are held constant
# thereafter. Users can of course create variable boundary conditions.
ij  = (slice(im), slice(jm))
ij0 = (slice(im), slice(jm), 0)
ua [ij]  = uab[ij]
va [ij]  = vab[ij]
el [ij]  = elb[ij]
et [ij]  = etb[ij]
etf[ij]  = et[ij]
d  [ij]  = h[ij] + el[ij]
dt [ij]  = h[ij] + et[ij]
w  [ij0] = vfluxf[ij]

ijk = (slice(im), slice(jm), slice(kb))
ij_ = (slice(im), slice(jm), None)
l   [ijk] = 0.1 * dt[ij_]   ##TODO Magic number!
q2b [ijk] = small
q2lb[ijk] = l[ijk] * q2b[ijk]
kh  [ijk] = l[ijk] * sqrt(q2b[ijk])
km  [ijk] = kh[ijk]
kq  [ijk] = kh[ijk]
aam [ijk] = aam_init

ijk = (slice(im), slice(jm), slice(kb-1))
q2 [ijk] = q2b [ijk]
q2l[ijk] = q2lb[ijk]
t  [ijk] = tb  [ijk]
s  [ijk] = sb  [ijk]
u  [ijk] = ub  [ijk]
v  [ijk] = vb  [ijk]

dens(s, t, rho)
baropg()

ij  = (slice(im), slice(jm))
k   = slice(kb-1)
ijk = (slice(im), slice(jm), slice(kb-1))
drx2d[ij] += dot(drhox[ijk], dz[k])
dry2d[ij] += dot(drhoy[ijk], dz[k])

# Calculate bottom friction coefficient
ij = (slice(im), slice(jm))
cbc[ij] = (kappa / log((1 + zz[-2]) * h[ij] / z0b)) ** 2
# If the following is invoked, then it is probable that the wrong choice
# of z0b or vertical spacing has been made:
cbc[ij] = clip(cbc[ij], cbcmin, cbcmax)

# Calculate external (2-D) CFL time step
ij = (slice(im), slice(jm))
tps[ij] = 0.5 / sqrt(1/dx[ij] ** 2 + 1/dy[ij] ** 2) / sqrt(grav * (h[ij] + small)) * fsm

# The following data are needed for a seamless restart. if nread=1,
# data had been created by a previous run (see write(71) at end of
# this program). nread=0 denotes a first time run.
##TODO Implement: lines 730--735

ij = (slice(im), slice(jm))
d [ij] = h[ij] + el[ij]
dt[ij] = h[ij] + et[ij]

time = time0

# Print geometry and other initial fields (select statements as desired)
##TODO Implement: lines 751--812

# Initial conditions. Select print statements in printall as desired
##TODO Implement: line 820

# Initialise netCDF output and output initial set of data
##TODO Implement: lines 826--829

# Begin internal (3-D) mode
for iint in xrange(iend):
    time = time0 + (iint + 1) * dti / spdi
    
    ##FIXME Incremental inertial ramp?
    if lramp:
        ramp = min(time / period, 1.0)
    else:
        ramp = 1.0
        
    # Set time dependent, surface and lateral boundary conditions.
    # The latter will be used in subroutine bcond. Users may
    # wish to create a subroutine to supply wusurf, wvsurf, wtsurf,
    # wssurf, swrad and vflux.
    
    ij  = (slice(1, im-1, 1), slice(1, jm-1, 1))
    ij0 = (slice(1, im-1, 1), slice(1, jm-1, 1), 0)
    
    e_atmos[ij] = 0.0
    vfluxf [ij] = 0.0
    
    # Set w(i,j,1)=vflux(i,j).ne.0 if one wishes non-zero flow across
    # the sea surface. See calculation of elf(i,j) below and subroutines
    # vertvl, advt1 (or advt2). If w(1,j,1)=0, and, additionally, there
    # is no net flow across lateral boundaries, the basin volume will be
    # constant; if also vflux(i,j).ne.0, then, for example, the average
    # salinity will change and, unrealistically, so will total salt.
    w[ij0] = vfluxf[ij]
    
    # Set wtsurf to the sensible heat, the latent heat (which involves
    # only the evaporative component of vflux) and the long wave
    # radiation.
    wtsurf[ij] = 0.0
    
    # Set swrad to the short wave radiation
    swrad[ij] = 0.0
    
    # To account for change in temperature of flow crossing the sea
    # surface (generally quite small compared to latent heat effect)
    ##FIXME Some sort of dark magic?
    tatm        = t[ij0] + tbias    # An approximation
    wtsurf[ij] += vfluxf[ij] * ((tatm - t[ij0]) - tbias)
    
    # Set the salinity of water vapor/precipitation which enters/leaves
    # the atmosphere (or e.g., an ice cover)
    satm       = 0.0
    wssurf[ij] = vfluxf[ij] * ((satm - s[ij0]) - sbias)
    
    